DALS008-数据探索(Exploratory Data Analysis)01-图表

发表于 | 分类于
| 字数统计: 9,044 | 阅读时长 ≈ 37

前言

这篇笔记是《Data Analysis for the Life Sciences》的第3章：探索性数据分析，主要内容涉及常见图表的解读，数据的描述，秩和检验。

图片的最大价值就在于能够强迫我们发现一些我们预料不到的东西——John W.Tukey

生命科学的研究，由于检测存在的误差，会导致一些异常值的出现。如果不没有发现这些问题，那么就会导致错误的分析结果以及假的发现。例如，有时候实验失败，或者是说在分析流程中，对涉及到的所有实验数据都进行检验，也无法发现它们，例如使用R中的t.test()函数。但是，这些分析流程仍然会给我们一些启示。此外，仅仅通过实验报告中的结果，也很难或者是几乎无法发现一些错误。

图表就是一种发现这些问题的强大工具，这一章节的内容就是探索性数据分析(Exploratory Data Analysis,EDA)。现存很多数据分析方法都是通过EDA来得到启示的。此外，EDA还能发现一些有意思的物理物现象，如果使用其它常规的统计学方法，这些现象有可能就会错失。在这本书中，我们会使用探索性数据来选择合适的统计学方法。

QQ图

qq图的全称为quantile quantile Plots，它经常用于验证一批数据是否服从正态分布。我们通过一些具体的案例来说明一下什么是分位数。

在一批数据中，第p个百分位数被定义为这个数大于p%的数，具体的定义可以参考以前的笔记《StatQuest学习笔记06——分位数及其应用》，再比如中间的50%百分位数就是中位数。

我们可以计算出一批身高的百分位数，如下所示：

1
2
3
4
library(UsingR) ##available from CRAN
library(rafalib)
x <- father.son$fheight
head(x)

结果如下所示：

1
2
> head(x)
[1] 65.04851 63.25094 64.95532 65.75250 61.13723 63.02254

现在我们计算一批服合正态分布的数据：

1
2
3
4
5
6
7
8
9
10
11
12
13
ps <- (seq(0,99) + 0.5 )/100
# Calculate standard quantile
qs <- quantile(x, ps)
# produces sample quantiles corresponding to the given probabilitie
# For example
# ps2 <- c(0.005,0.515,0.995)
# quantile(x,ps2)
# > quantile(x,ps2)
# 0.5%    51.5%    99.5% 
# 60.29272 67.84169 74.58332 
normalqs <- qnorm(ps, mean(x), popsd(x))
plot(normalqs,qs,xlab="Normal percentiles",ylab="Height percentiles")
abline(0,1) ##identity line

上图显示的是，(mean(x))为均值，总体方差为（pposd(x)）的正态分布数据集的qq图，这里的x指的是前文提到的father.son中fheight的身高，也就是父亲的身高，这里的qq图是理论上的qq图，代码中的qnorm()函数用于计算正态分布中的第p个分位一数的Z-score，具体的可以参考这篇笔记《R语言中dnorm, pnorm, qnorm与rnorm以及随机数》，现在看一下实际的qq图，如下所示：

1
2
qqnorm(x)
qqline(x)

qqnorm()函数绘制的图形都是按照标准正态分布来进行绘制，因此，上图中的直接斜率为popas(x)，截矩为mean(x)。

在前面的案例中，多数点与直线的匹配程度都非常好，实际上，我们可以使用Monte Carlo来模拟一下已知数据是否符合正态分布，如下所示：

1
2
3
4
n <-1000
x <- rnorm(n)
qqnorm(x)
qqline(x)

现在我们来模拟一下不符合正态分布的数据在qq图中是什么样子的，我们从不同自由度的t分布的数据中来抽取一些数据来画图，服从t分布的数据有这么一个特点：自由度越低，t分布的尾部就越平坦，自由度越高，尾部越高，整体越接近正态分布。在qq图中，我们可以看到，在qq图的左侧，这些数据点低于直线，而在右侧，这些点高于直线，如下所示：

1
2
3
4
5
6
7
dfs <- c(3,6,12,30)
mypar(2,2)
for(df in dfs){
x <- rt(1000,df)
qqnorm(x,xlab="t quantiles",main=paste0("d.f=",df),ylim=c(-6,6))
qqline(x)
}

箱线图(boxplot)

我们获取的数据并非总是符合正态分布的，例如收入(income)就不符合正态分布，对于那些不符合正态分布的数据来说，谈论均值与标准差就没有多大意义，因为我们无法从均值与标准差中推测出相应的总体信息，只有数据符合正态时，谈论均值与标准差才有意义，现在我们来看一下2000年美国CEO的收入，如下所示：

1
2
3
4
5
6
library(rafalib)
library(UsingR)
mypar(1,2)
hist(exec.pay)
qqnorm(exec.pay)
qqline(exec.pay)

除了使用QQ外，另外一种估计数据的实用手段就是计算一组数据的25%，50%与75%分位数。箱线图就能展示这3个计算值，以及这些数据的范围，即中位数±1.5（75%百分位数-25%百分位数）。超过这个范围的数据点有时候称为异常值(outliers)。

1
boxplot(exec.pay, ylab="10,000s of dollars",ylim=c(0,400))

散点图与相关(Scatterplots And Correlation)

前面提到的方法都是单变量(univariate variables)。而在生物学研究中，我们往往关注两个变量或多个变量之间的关系。一个典型的案例就是Francis Galton通过计算父亲与儿子身高的关系来理解遗传学。如果我们想要计算一下这些数据，我们就能利用这两组数据（也就是父亲的身高这一组数据与儿子的身高这一组数据）的均值与标准差来描述这些数据，因为身高这种数据是符合正态分布的。但是，我们通常这种常规方法无法找出这两组数据的关系，先来看一下散点图，如下所示：

1
2
3
4
5
library(UsingR)
data("father.son")
x <- father.son$fheight
y <- father.son$sheight
plot(x,y, xlab="Father's height in inches",ylab="Son's height in inches",main=paste("Correlation=", signif(cor(x,y),2)))

这个散点图展示出了一个趋势：父亲越高，儿子就越高。这个趋势可以使用相关系数(correlation coefficient)来表示，这个案例中的相关系数是0.5。也就是说我们可以使用父亲的身高来预测儿子的身高。

分层(Stratification)

现在我们假设一种场景，比如我们随机地选取一个儿子来猜测他的身高，其中儿子这一组数据的平均身高是68.7英寸（），这个身高所占的比例是最大的（可以看直方图），因此我们可以根据这个预测一个儿子的身高。但是，如果我们被告知，这个选中的儿子的父亲身高是72英寸，那么我们还会猜测这个儿子身高是68.7英寸么？

父亲的身高是高于平均值的，尤其是当这个父亲的身高还高于父亲这一组数据均值的1.75个标准差时。我们有可能会预测，这个父亲的儿子身高是否也高于1.75个标准差？结果发现，这个结果是高估了。为了验证这个结果，我们来看一下，那些身高约为72英寸的父亲的儿子的身高，我们可以对父亲这一组数据进行分层(stratifying)。

1
2
3
4
groups <- split(y, round(x))
boxplot(groups)
print(mean(y[ round(x) == 72]))
# [1] 70.67719

对数据进行分层后，再绘制箱线图，我们就可以看到每组数据的分布。身高约为72英寸父亲的平均身高是70.7英寸。我们还能看到，分层后的数据的中位数大概是一条直接（箱线图中间的线是中位数值，而非均值）。中位数构成的这条直线比较类似于回归线(regression line)，而回归线的斜率与相关有关。

二元正态分布(Bi-variate Normal Distribution)

一对随机变量，例如(X,Y)，当它们值的比例可以按照以下公式来表示时，例如a和b，就可以认为它们接近二元正态分布：

上面的这个公式看起来很复杂，但是背后的原理非常简单。另外一种比较简单的理解就是：先固定一个值x，然后观察当X=x时，所有的点(X,Y)。通常来说，在统计学中我们把这一过程称为设定条件(conditioning)。我们可以根据X来设定Y。如果一对随机变量接近双变量正态分布，那么X=x时，Y的分布也就服从一个正态分布，而x的值可以是任意的，现在我们来看一下身高数据，看一下上述的思路是否符合，身高数据一共分为4层，如下所示：

1
2
3
4
5
6
7
groups <- split(y, round(x))
mypar(2,2)
for (i in c(5,8,11,14)){
  qqnorm(groups[[i]], main = paste0("X=", names(groups)[i]," strate"),
         ylim=range(y), xlim=c(-2.5, 2.5))
  qqline(groups[[i]])
}

当两个变量服从双变量正态分布时，从数学统计学的公式就可以知道，对于给定的x值，相对应的Y（与X=x配对）平均值为：

从上面公式我们可以知道，这个直线的斜率为$\rho \frac{\sigma_{Y}}{\sigma_{X}}$，这个直线就是回归直线(regression line)。如果SD一样，那么这个回归线的斜率就是$\rho$。因此，如果我们对X和Y进行标准化后，那么相关(correlation)就是这个回归直线的斜率。

这个公式的另外一种表示就是我们可以构建一个预测值$\hat{Y}$：即对于距离x均值1个SD的x值，我们可以预测相应Y的值，是远离y均值的$\rho$个SD的值，如下所示：

如果两个变量，即X和Y存在着很好的相关性，那么我们可以预测出相同的SD数目。如果相关性为0，那么我们就无法根据x的值来预测y的值。现在我们在计算一下每层的均值，进而确定回归直线：

1
2
3
4
5
6
7
x=( x-mean(x) )/sd(x)
y=( y-mean(y) )/sd(y)
means=tapply(y, round(x*4)/4, mean)
fatherheights=as.numeric(names(means))
mypar(1,1)
plot(fatherheights, means, ylab="average of strata of son heights", ylim=range(fatherheights))
abline(0, cor(x,y))

方差解释(Variance explained)

前面提到的条件分布(conditional distribution)的标准差如下所示：

这个公式与X能够解释Y的多少方差类似，这个解释程度是$\rho^{2} \times 100 \%$。而上述的那个公式中：Y的方差是$\rho^{2}$，如果给定X的值，那么Y自身能够解释Y值波动的程度就是$
\left(1-\rho^{2}\right) \sigma_{Y}^{2}
$。这个里我们只需要记住，方差解释(variance explained)这种表述只有在两组数据接近双变量正态分布(bivariate normal distribution)的情况下才有意义（注：对于书中的这段话我没太明白，按照自己的理解来记的笔记）。

绘图原则

优秀数据图主要是为了能够准确清楚地展示你的数据，Karl Broman（注：此人是华盛顿大学生物统计与医学信息学系教授，主要研究统计遗传学）提出了几点绘图中要避开以下操作：

  • 展示信息太少(Display as little information as possible)；
  • 阅读体验不好（乱堆图表）(Obscure what you do show (with chart junk))；
  • 使用伪3D图形，瞎用颜色(Use pseudo-3D and color gratuitously)；
  • 随意使用饼图（尤其是各种颜色与3D饼图）(Make a pie chart (preferably in color and 3D))；
  • 没有图例(Use a poorly chosen scale)；
  • 没有展示有效数字(Ignore significant figures)。

常见各种图形

这一部分涉及常见的各种图形，例如饼图，柱状图等。

饼图(pie charts)

现在我们来看一个数据，这个数据是2013年的调查的各个浏览器的使用率，如下所示：

1
2
browsers <- c(Opera=1,Safari=9,Firefox=20,IE=26,Chrome=44)
pie(browsers,main="Browser Usage (August 2013)")

但是饼图有很大局限，在R中，我们查看一下pie()函数的说明，如下所示：

“Pie charts are a very bad way of displaying information. The eye is good at judging
linear measures and bad at judging relative areas. A bar chart or dot chart is a preferable
way of displaying this type of data.”

上述就提到了，饼图展示的信息有限。人眼善于判断线性度量，不善于判断面积度量。因此，使用条形图来展示相关的信息比饼图更好。

从饼图中我们只能看出每个浏览器的占有率，但是，除非这些占有率接近25%，50%或75%，否则也不太容易看出来，此时如果直接展示数字反而好一些。为了更好地展示数据，现在使用柱状图(barplot)来看一下这些数据，如下所示：

1
2
3
barplot(browsers, main="Browser Usage (August 2013)", ylim=c(0,55))
abline(h=1:5 * 10)
barplot(browsers, add=TRUE)

从柱状图上的柱子，以及平行于x轴的几条水平线我们就很容易看到每个浏览器的市场占有率。不过在画柱状图的时候，要避免画成3D的，因为3D图很直观地看到其顶部对应的数字，如下所示：

更差的是画环形图(donut plot)，如下所示：

柱状图(barplot)

这里注明一下，barplot有的翻译为柱状图，有的翻译为条形图，我觉得是一样的，这两种图形的x轴（或y轴）上都是分类变量，而不是连续变量。需要与柱状图进行区分的则是直方图(histogram)，直方图的x轴上是连续变量。

虽然柱状图在展示百分比方面很有用，但是柱状图经常被错误地用于比较两组数据，在这类应用中，常常用均值作为柱子的高度，使用误差棒来表示标准误，如下所示：

这种图展示的信息很少，例如没有说明数据的分布情况。其实针对这种情况，使用箱线图(boxplot)更合适，如果数据点比较少，可以直接在箱线图上添加上数据点。如果数据点很多，只用画箱线图，不加数据点注行。如果在箱线图函数boxplot()中添加上range=0参数的话，也可以不绘制离群点(outlier)，现在我们绘制箱线图，如下所示：

1
2
3
4
5
6
7
8
x <- rnorm(8,30,5)
y <- rnorm(8,38,5)
library(rafalib)
mypar()
dat <- list(Treatment=x,Control=y)
boxplot(dat,xlab="Group",ylab="Response",cex=0)
stripchart(dat,vertical=TRUE,method="jitter",pch=16,add=TRUE,col=1)

从上面的箱线图我们可以看到以下信息：

  • 数据的中心；
  • 数据的分布，范围；
  • 所有的数据点。

而在柱状图中，我们仅能看到均值与SE，况且SE还与样本的数目(sample size)有关，无法发现数据的分布规律。当我们的数据有离群点或者是有很大拖尾时(very large tails)时，柱状图的一些局限就会更加明显，在下面的图形中我们会觉得这两组数据的差异非常大：

但是，如果我们仔细研究一下这两组数据，就会发现，这两组数据的差异主要是由第二数据的两个离群值导致的，为了解决这个问题，现在我们分别绘制一下常规的箱线图以及经过log转换后的箱线图，如下所示：

1
2
3
4
5
6
7
8
9
10
11
library(rafalib)
mypar(1,2)
x <- 10^rnorm(10,0.2)
y <- 10^rnorm(10,1.2)
dat <- list(Treatment=x,Control=y)
boxplot(dat,xlab="Group",ylab="Response",cex=0)
stripchart(dat,vertical=TRUE,method="jitter",pch=16,add=TRUE,col=1)
boxplot(dat,xlab="Group",ylab="Response",log="y",cex=0)
stripchart(dat,vertical=TRUE,method="jitter",pch=16,add=TRUE,col=1)

散点图(scatter plot)

许多统计分析研究的目标就是找到两个变量的关系。在这些统计结果中通常会报告简单的相关性(correlation)，有时候还会画出相应的图形。但是，仅仅展示出一个回归直线(regression line)还不够，因为直线无法显示出背后的数据，此时最好还要加上数据点，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
x <- rnorm(100,100,10)
y <- 100+10*sqrt(2*(x-min(x)))+rnorm(100,0,10)
mypar(1,2)
plot(x,y,lwd=2,type="n")
fit <- lm(y~x)
abline(fit$coef,lwd=2)
b <- round(fit$coef,4)
text(78, 200, paste("y =", b[1], "+", b[2], "x"), adj=c(0,0.5))
rho <- round(cor(x,y),4)
text(78, 187,expression(paste(rho," = 0.8567")),adj=c(0,0.5))
plot(x,y,lwd=2)
fit <- lm(y~x)
abline(fit$coef,lwd=2)

高度的相关性并不意味着相同(High correlation does not imply replication)

当出现一个新的技术，或者是在实验室中引入一项新技术时，我们通常会用这些技术来检测相同的样本(replicated sample)，进而绘制出散点图来看一下它们的相关性。高度的相关性(high correlation)用于评估这项新技术是否有重现性(reproducible)。

但是，相关性有时候也会误导人(misleading)，下面的散点图展示的是一项高通量技术对相同样本的检测结果。这项技术同时输出12626个检测值。在下图的左侧中，我们可以看到原始数据显示出很高的相关性。但是这组数据拖尾严重(very fat tails)，其中有95%的数据是在绿线以下，另外，下图的右侧是经过log转换后的图：

虽然上面两张图都显示出了高度的相关性（都接近1）。这是否就意味着这些数据能够重复？为了研究第二次检测的结果能否重现第一次的检测，我们需要研究这两组数据的差值。为此，我们可以绘制差值与均值的图型，如下所示：

上面的这种分析方法叫Bland-Altman绘图法，也叫MA图。MA的意思是minus与average，因为在这种图中上，y轴表示了两个样本数据经过log转换后，取的差值，x轴表示的是两个样本的数据经过log转换后，两个数据的平均值。在这张图形中，我们使用的是log2转换，从y轴上可以看到，两次测量后的差值多数情况下是1。这意味着当测量应该相同时，我们将平均观察到2倍的差异。我们现在可以将这种可变性(variability)与检测的差值(differences)进行比较，从而决定这项技术是否达到了我们想要的精度。

配对数据的柱状图

数据分析的一项常规任务就是比较两组数据。当数据集比较小，并且是配对数据（例如治疗前与治疗后），可以使用两种颜色的柱状图来展示结果，如下所示：

除此之外，还有更好的展示数据的方式，这些展示方式能够看到治疗效果。第一种就是简单地绘制散点图，从散点图中我们可以发现，多数点都在直线上方。另外一种就是绘制箱线图。
先看散点图，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
set.seed(12201970)
before <- runif(6, 5, 8)
after <- rnorm(6, before*1.05, 2)
li <- range(c(before, after))
ymx <- max(abs(after-before))
mypar(1,2)
plot(before, after, xlab="Before", ylab="After",
     ylim=li, xlim=li)
abline(0,1, lty=2, col=1)
plot(before, after-before, xlab="Before", ylim=c(-ymx, ymx),
     ylab="Change (After - Before)", lwd=2)
abline(h=0, lty=2, col=1)

从上图可以看出来，右侧的图比较像MA图，其实就是相当于把左侧的散点图旋转了一下。

散点图中的直线并不直观，而箱线图则能很好地显示变化的效果，但是会丧失配对信息，如下所示：

1
2
3
4
5
6
7
8
z <- rep(c(0,1), rep(6,2))
mypar(1,2)
plot(z, c(before, after),
xaxt="n", ylab="Response",
xlab="", xlim=c(-0.5, 1.5))
axis(side=1, at=c(0,1), c("Before","After"))
segments(rep(0,6), before, rep(1,6), after, col=1)
boxplot(before,after,names=c("Before","After"),ylab="Response")

折线图(line chart)

先看一个反面典型：

这个折线图是一个伪3D图，它表示的信息不清楚，也许是为了区分这三条曲线，但是，从这3条曲线上的点很难找到对应的数值，这种图不推荐使用。

现在来看一个比较正常的图表，代码如下所示：

1
2
3
4
5
6
7
8
9
10
log_dose <- 1:10
Drug_A <- c(0.8, 0.75, 0.76, 0.71, 0.73, 0.7,0.67,0.66,0.68,0.65)
Drug_B <- c(0.75,0.73,0.7,0.65,0.62,0.5,0.3,0.28,0.27,0.29)
Drug_C <- c(0.92,0.9,0.87,0.75,0.74,0.73,0.64,0.6,0.58,0.52)
x <- data.frame(log_dose,Drug_A,Drug_B,Drug_C)
plot(x[,1],x[,2],xlab="log Dose",ylab="Proportion survived",ylim=c(0,1),
     type="l",lwd=2,col=1)
lines(x[,1],x[,3],lwd=2,col=2)
lines(x[,1],x[,4],lwd=2,col=3)
legend(1,0.3,c("Drug A","Drug B","Drug C"),lwd=2, col=1:3,horiz=T,bty = "n")

添加重要信息

在下面的案例中，我们会模拟生成一批数据来研究两组之间的关系：对照组与治疗组。每组分别有3次实验，现在我们使用普通的柱状图来描述一下这两种干预结果，如下所示：

这类图形并没有很好地展示不同药物剂量的治疗效果，也就是没有展示出量效关系这种重要信息，现在我们使用折线图来展示一下这些信息，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
x <- 1:8
y1 <- 0.9 - x/80 + rnorm(length(x), 0, 0.02)
y2 <- 0.9 - x/40 + rnorm(length(x), 0, 0.02)
y3 <- 0.85 - x/30 + rnorm(length(x), 0, 0.02)
y <- cbind(y1, y2, y3)
z1 <- 0.95 - x/40 + rnorm(length(x), 0, 0.02)
z2 <- ilogit(-0.4*(x-4.5) + rnorm(length(x), 0, 0.04))
z3 <- ilogit(-0.5*(x-4.5) + rnorm(length(x), 0, 0.04))
z1[6:8] <- z1[6:8] - 0.18*3
z <- cbind(z1, z2, z3)
ym <- apply(y, 1, mean)
zm <- apply(z, 1, mean)
dev.off()
png(file="../fig9b.png", width=1200, height=1000, res=288,pointsize=12)
bgcolor <- rgb(0, 0, 98, maxColorValue=255)
par(las=1,fg="white",col="white",col.axis="white",col.lab="white",
    bg=bgcolor,mar=c(0.1,0.1,0.1,0.1))
plot(x, y1, ylim=c(0,1), type="n", xlab="Dose", ylab="Response")
for(i in 1:3) lines(x, y[,i], col=1, lwd=1, lty=2)
for(i in 1:3) lines(x, z[,i], col=2, lwd=1, lty=2)
lines(x, ym, col=1, lwd=2)
lines(x, zm, col=2, lwd=2)
legend("bottomleft", lwd=2, col=c(1, 2), c("Control", "Treated"))
dev.off()

有效数字

R的统计结果的小数点数位很多，有的时候直接粘贴比较麻烦，此时就可以减少一些小数位，先来看一批数据，如下所示：

1
2
3
4
5
heights <- cbind(rnorm(8,73,3),rnorm(8,73,3),rnorm(8,80,3),
                 rnorm(8,78,3),rnorm(8,78,3))
colnames(heights)<-c("SG","PG","C","PF","SF")
rownames(heights)<- paste("team",1:8)
heights

结果如下所示：

1
2
3
4
5
6
7
8
9
10
> heights
             SG       PG        C       PF       SF
team 1 70.28569 73.53063 77.97624 75.53098 79.91438
team 2 69.13691 72.98393 79.53870 77.05309 73.23242
team 3 69.66017 70.11667 80.43766 80.39365 76.69519
team 4 73.95542 72.61165 77.36032 79.48383 75.00610
team 5 73.36905 77.19182 78.59511 78.28309 81.98515
team 6 74.03342 67.47227 79.25801 70.24398 77.07891
team 7 74.14659 76.21967 81.25641 74.01896 80.90877
team 8 70.57122 68.65044 84.30345 76.18054 77.83099

上述的这些数字精确到小点第5位，很多时候在写实验结果时根本没必要，此时使用round()函数，如下所示：

1
2
3
4
5
6
7
8
9
10
> round(heights,1)
         SG   PG    C   PF   SF
team 1 70.3 73.5 78.0 75.5 79.9
team 2 69.1 73.0 79.5 77.1 73.2
team 3 69.7 70.1 80.4 80.4 76.7
team 4 74.0 72.6 77.4 79.5 75.0
team 5 73.4 77.2 78.6 78.3 82.0
team 6 74.0 67.5 79.3 70.2 77.1
team 7 74.1 76.2 81.3 74.0 80.9
team 8 70.6 68.7 84.3 76.2 77.8

相关的误解(Misunderstanding Correlation)

相关性(correlation)广泛用于表示数据的重现性(reproducibility)，例如在基因组学(genomics)常会用到。但在数学角度来看，相关性与重现性关系不大，这里简单描述三个关于相关性的问题。

关于相关性最大的误用就是计算那些不服从双变量正态分布数据的相关性。前面提到，如果两组数据服从双变量正态分布，那么使用均值、标准差以及相关性来描述这些数据是合适的。但是，还有很多数据不是双变量正态分布。例如基因表达的数据就不服从双变量正态分布，这种数据的右侧有着扁平的尾巴(with a very fat right tail)。

研究两个重复检测数据重现性(reproducibility)的标准做法是先计算一下这两组数据的距离，如下所示：

当两次检测的数值重复性很好的时候，这个度量(metric)是降低的，当这个值是0的时候（其实就相当于两次检测的结果完全一样），结果最完美（这只是理想情况）。使用这种度量的另外一个好处就是，如果我们将它除以N，那么得到的数值就是就是作为$d_{1},\dots,d_{N}$的标准差（此时我们假设d是大于0的）（原文是if we assume the d average out to 0，我的理解是，我们假设d的平均值为0，或者是d的平均值大于0）。如果将d视为残差(residuals)，那么这个统计量（已经除以N了）就相当于RMSE（全称为Root Mean Squared Error，译为均方根误差），我们称为距离度量(distance metric)。此外，这个统计量的单位与检测值的单位是一样的，因此能很好地解释结果。

使用相关性(correlation)是另外一个局限是，相关性计算无法检测到由于平均值改变而导致的不可重复性（原文是：Another limitation of the correlation is that it does not detect cases that are not reproducible due to average changes.），而距离度量则能检测到这些差异，因此前面的公式可以写为下面的形式：

其中$\mu_{x}$和$\mu_{y}$表示的是两次检测的均值，因此上面的公式还可以写为如下形式：

如果我们假设前后两次检测的结果的方差都是1，那么上面的公式就可以简化为如下形式：

其中$\rho$表示相关性，因此我们可以直接看到距离度量(distance metric)与相关性(correlation)的关系，但是，距离度量与相关性的一个重要差异就是距离度量含有$\left(\mu_{x}-\mu_{y}\right)^{2}$，因此，它可以检测到那些由于大量均值变化而导致的无法重复。

相关性不能当成重复性(reproducibility)指标的另外一个因素就是，相关性没有单位(units)。我们使用一个公式就能看到这种局限，假如我们有两个变量，分别为x和y，其中y=x+d。如果d的方差越大，那么x+d就越无法重现x，也就是说，如果d这个数值波动越大，表示的就是第二次检测与第一次检测的数值偏差越大，进而可以得出结果：两次检测数据之间的重复性越差。而距离度量(distrance metric)仅取决于d的方差，也就是说它只取决于前后两次检测的偏差程度，它表示的是就是重复性的好坏，而相关性还依赖于x的方差，如果d与x是独立的，那么相关性的公式就如下所示：

从上面的公式我们就可以看出来，相关性取决于d的方差与x的方差，这就说明了，即使相关性接近于1，也无法说明重复性很好。具体来说就是，无论d的方差怎么样，我们通过增加x的方差，就可以使相关系数随意地接近于1。

练习

P131

异常值的总结

正态分布在分析生命科学相关数据中有着重要作用，但是，由于相关检测仪器的复杂性，也会偶尔观察到一些异常值。例如某些扫描仪的缺陷会出现几个异常值的点，或者说是PCR过程中会产生一些偏倚。

因此我们往往会遇到一些情况，即在一群数据中，除了某个异常值外，其余的数据符合正态分布，如下所示：

1
2
3
4
set.seed(1)
x <- c(rnorm(100,0,1))
x[23] <- 100
boxplot(x)

例如在上面的图形中，有一个点的数值是100，在统计学中，我们将这种点称为异常值(outliers)。一小部分的异常值会严重误导我们的数据分析，例如，这组数据中的均值与方差差就会远远偏离实际值，如下所示：

1
2
> cat("The average is",mean(x),"and the SD is",sd(x))
The average is 1.108142 and the SD is 10.02938

我们生成的数据是标准正态分布，也就是均值为1，标准差为0，而加了一个异常值后，均值就变成了1.108142，标准差就变成了10.02938。

中位数

中位数的定义就是，将一组数据按从小到大的顺序排列，这个数字位于中间，有50%的数字大于它，有50%的数字小于它，对于含有异常值的一组数据来说，中位数比较稳健(robust)，现在我们来看一下我们前面数据的中位数（理论上应该比较接近于0），如下所示：

1
median(x)

结果如下所示：

1
2
> median(x)
[1] 0.1684483

中位数绝对偏差(The median absolute deviation,MAD)

中位数绝对偏差，即The median absolute deviation，缩写为MAD，它是对标准差的稳健性描述。它的计算方法为：

第一步：计算每个数据点与中位数的差值；

第二步：取第一步中的绝对值；

第三步：取第二步中绝对值的中位数，再乘以1.4826。

其中，1.4826是一个比例因子(scaling factor )，MAD是一个无偏估计，现在我们使用MAD来估计前面的那组数据（就是服从标准正态分布的那组数据），那么它的值就非常接近于1了，如下所示：

1
2
> mad(x)
[1] 0.8857141

Spearman相关系数

前面我们提到，相关性对异常值比较敏感，这里我们构建了两组独立的数据，基中它们都含有以前面类似的异常值，如下所示：

1
2
3
4
5
6
7
8
9
set.seed(1)
x=c(rnorm(100,0,1)) ##real distribution
x[23] <- 100 ##mistake made in 23th measurement
y=c(rnorm(100,0,1)) ##real distribution
y[23] <- 84 ##similar mistake made in 23th measurement
library(rafalib)
mypar()
plot(x,y,main=paste0("correlation=",round(cor(x,y),3)),pch=21,bg=1,xlim=c(-3,100),ylim=c(-3,100))
abline(0,1)

Spearman相关系数遵循中位数与MAD的基本思想，也就是使用了分位数(quantiles)的原理。这个原理很简单：我们将每个数据集转换为秩(rank)，然后计算它们的相关性，如下所示：

1
2
3
4
mypar(1,2)
plot(x,y,main=paste0("correlation=",round(cor(x,y),3)),pch=21,bg=1,xlim=c(-3,100),ylim=c(-3,100))
plot(rank(x),rank(y),main=paste0("correlation=",round(cor(x,y,method="spearman"),3)),pch=21,bg=1,xlim=c(-3,100),ylim=c(-3,100))
abline(0,1)

从上面的结果我们可以知道，这些统计量对异常值很稳健，但为什么我们会使用那些不稳健的相关性计算方法呢？通常来说，如果我们知道异常值，那么我们最好使用中位数和MAD来替代均值与标准差来描述一组数据。但是，与非稳健性的相关性计算方法相比，这些稳健的统计方法的功效比较低。

对数比的对称性(Symmetry of log ratio)

比值(ratio)通常并不对称，现在我们来模拟一组比值来看一下，这组比值表示的是两个样本中100个基因的比值，如下所示：

1
2
3
x <- 2^(rnorm(100))
y <- 2^(rnorm(100))
ratios <- x / y

在生命科学中，我们经常使用比值或表达倍数(fold changes, FC)。假设现在们就在研究某个因素干预前后基因的变化的比值。因此在这些数据中，有些比值大于1，就表示经过干预后，基因表达上升。如果干预没有效果，那么这些比值小于1和大于1的数目是差不多的，通过直方图我们就可以看到干预后的基因变化量，如下所示：

1
2
3
4
mypar(1,2)
hist(ratios)
logratios <- log2(ratios)
hist(logratios)

从上面的这两张图我们可以看到，这些数据并不关于1对称。例如，1/32比32/1更接近于1，其实这两个数字表示的比值意义其实是相反的，一个是降低了32倍，一个是升高了32倍，但是它们不对称。如果我们经过了log转换，这两个数值就会关于0对称，转换的过程如下所示：

现在我们看一下这些数据点，如下所示：

在生命科学中，经常使用log转换，这是因为100倍的变化可能是100/1或1/100，但是1/100比100/1更接近于1。

Wilcoxon秩和检测(Wilcoxon Rank Sum Test)

从前面内容我们可以知道，样本的均值(mean)和SD对异常值(outlier)敏感。t检验的计算原理就是均值与SD。因此在有异常值存在的情况下，Wilcoxon秩和检验（其本质与Mann-Whitney检验）就成了一种替代选择。

在下面的代码中，我们进行了t检验（其零假设为真）。当使用这些原始数据进行计算时，t检验与Wilcoxon检验的p值都大于0.05，但是，当我们在每个样本中故意错误地更改了一个观测值后，并将其输入计算过程，此时，t检验就会计算出一个小p值，而Wilcoxon检验则不会这样，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
set.seed(779) 
##779 picked for illustration purposes
N=25
x<- rnorm(N,0,1)
y<- rnorm(N,0,1)
cat("t-test pval:",t.test(x,y)$p.value,"\n")
cat("Wilcox test pval:",wilcox.test(x,y)$p.value,"\n")
x[1] <- 5
x[2] <- 7
cat("t-test pval:",t.test(x,y)$p.value,"\n")
cat("Wilcox test pval:",wilcox.test(x,y)$p.value,"\n")

计算结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
> cat("t-test pval:",t.test(x,y)$p.value,"\n")
t-test pval: 0.1589672 
> cat("Wilcox test pval:",wilcox.test(x,y)$p.value,"\n")
Wilcox test pval: 0.3066718 
> x[1] <- 5
> x[2] <- 7
> cat("t-test pval:",t.test(x,y)$p.value,"\n")
t-test pval: 0.04439948 
> 
> cat("Wilcox test pval:",wilcox.test(x,y)$p.value,"\n")
Wilcox test pval: 0.1310212

Wilcoxon检验的基本思想

Wilcoxon检验的基本思想如下所示：

  1. 汇总所有数据；
  2. 将数据的值转换为秩(rank)，也就是排序；
  3. 然后将这些秩再按原来的组别分开；
  4. 计算出每组秩的总和，进行检验。

检验的过程如下所示：

1
2
3
4
5
6
7
8
9
library(rafalib)
mypar(1,2)
stripchart(list(x,y),vertical=TRUE,ylim=c(-7,7),ylab="Observations",pch=21,bg=1)
abline(h=0)
xrank<-rank(c(x,y))[seq(along=x)]
yrank<-rank(c(x,y))[-seq(along=y)]
stripchart(list(xrank,yrank),vertical=TRUE,ylab="Ranks",pch=21,bg=1,cex=1.25)
ws <- sapply(x,function(z) rank(c(z,y))[1]-1)
text( rep(1.05,length(ws)), xrank, ws, cex=0.8)

上图是含有异常值的两个图表，左侧是原始图形，它含有的异常值，右则是将数值点转换为秩后的图形，可以发现右侧的数据分布就比较正常了，没有异常值。

秩和检验中的W值是对第一组相对于第二组的秩。通过排列组合(combinatorics)，我们可以精确地计算出W的p值。我们也可以使用CLT来计算，因为W是近拉于正态分布的，下面我们构建了z-score来进行计算，如下所示；

1
2
3
4
5
W <-sum(ws)
W
n1<-length(x);n2<-length(y)
Z <- (mean(ws)-n2/2)/ sqrt(n2*(n1+n2+1)/12/n1)
print(Z)

结果如下所示：

1
2
3
4
5
6
7
> W <-sum(ws)
> W
[1] 391
> n1<-length(x);n2<-length(y)
> Z <- (mean(ws)-n2/2)/ sqrt(n2*(n1+n2+1)/12/n1)
> print(Z)
[1] 1.523124

在这个计算过程中，我们得到了Z值，这个Z值对应的p值大于0.05。

练习

P140